\documentclass[aps,pre,a4paper]{revtex4}
\usepackage{geometry}                % See geometry.pdf to learn the layout options. There are lots.
\usepackage{graphicx}
\usepackage{amssymb}
\usepackage{epstopdf}
\usepackage[utf8]{inputenc} % pour les accents sur les macs
%\usepackage[latin1]{inputenc} % pour les accents sur les pcs.
\usepackage[english]{babel}
\usepackage{amsmath}
\usepackage{color}
\usepackage{pstricks}
\usepackage[section]{placeins}
\DeclareGraphicsRule{.tif}{png}{.png}{`convert #1 `dirname #1`/`basename #1 .tif`.png}

%\usepackage{pslatex}
%\usepackage{bookman}


\begin{document}
\title{Verlet and Beeman algorithms}
\author{Bruno Garbin}
\date{\today}

\maketitle
\tableofcontents %table des matières
\newpage % on change de pages


\section{Generalities}
We will describe here Velocity Verlet and Beeman algorithms and their way of using in molecular dynamics simulation.

In molecular dynamics simulation we would simulate the movement of particles in a given place, so several points are very important to obtain physics results, like: the interaction between particles, the environment of simulation, of course one algorithm which calculate the different part of movment (position, velocity, acceleration)...

To begin it's necessary to choose the conditions of simulation, like the size of the system that we would simulate or still the initial position of particles. 
We want simulate an infinite system, but we can't take an infinite number of particles. So we must take periodic bounding conditions (PBC) around one box, to "reproduce" an infinite system. In fact that amounts at consider an infinite box's number, as it's shown on the following figure. 
$\\ \\ \\ \\$
%\begin{figure}[htbp] %  figure placement: here, top, bottom, or page
\begin{center}
	   \includegraphics[width=0.5 \columnwidth]{figure1.pdf} 
	   \label{fig:figure1}
\end{center}
%\end{figure}
$\\ \\ \\ \\$

When the particle $\alpha$ arrive to the right bound, it is "teleported" to the otherside of the box and go from position \fbox{1} to the position \fbox{2}.
This fact imply that the number of particles stays constant because a particle which go out the box enter by the otherside.

It's need to consider after the interaction between particles which imply the movment, we take here the Van Der Waals potential which is attractive if the two  particles considered are distant and repulsive if the two particles are too close (less than $\sigma$ the distance in the equilibrium). 
To this potential calcul it is necessary to consider each pair of particles that is to calculate the force applied at the particle $\alpha$ we must calculate the potential that the particle $\beta$ applied its, then the potential that the particle $\gamma$ applied its... for all particles.

%\begin{figure}[htbp] %  figure placement: here, top, bottom, or page
\begin{center}
	   \includegraphics[width=0.5 \columnwidth]{figure2.pdf} 
	   \label{fig:figure2}
\end{center}
%\end{figure}
$\\$
$\\$
$\\$
$\\$
$\\$
Don't forget the PBC,by exemple if the particle $\alpha$ is closer to the particle $\delta$ in the box 3 than if the particle $\delta$ was in the real box, we must consider the potential associate a the case 3. So to calculate the potential we consider all particles in the red box (which is the same dimension than real box but centered in $\alpha$).
$\\ \\ \\ \\$

%\begin{figure}[htbp] %  figure placement: here, top, bottom, or page
\begin{center}
	   \includegraphics[width=0.5 \columnwidth]{figure3.pdf} 
	   \label{fig:figure3}
\end{center}
%\end{figure}
$\\ \\ \\$
\section{Potentials}
We owe considered one potential, which will calculate as explain in previous part. So we take two different potential for the moment, which are Lennard-Jones potential and Lennard-Jones potential truncated and shifted. We can observe its on the following figure, for $\epsilon$ and $\sigma$ equal in one: 

%\begin{figure}[htbp] %  figure placement: here, top, bottom, or page
\begin{center}
	   \includegraphics[width=0.5 \columnwidth]{graphe.jpg} 
	   \label{fig:figure4}
\end{center}
%\end{figure}

\subsection{Lennard-Jones potential}
The "normal" Lennard-Jones potential has for expression:

\begin{equation}
	\boxed{V_{LJ}(r) = 4 \epsilon \left( \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right)}
\end{equation}

It is attractive at long range (this term correspond at Van-Der-Waals forces) and repulsive at short range (what corresponds at Pauli repulsion). As we can see the potential beetween two particles very distant is not zero, so when we calculate the force which is applied at one, it's need to consider all. When we want simulate this, if we have lots of particles, even with a successful computer, it's \textbf{very} long to calculate forces .

\subsection{Truncated and shifted Lennard-Jonnes potential without force discontinuity}
The truncated and shifted Lennard-Jones potential has for expression:

\[
\left \{
\begin{array}{c @{=} c}                                 
	V_{LJ}^{ts}(r) \ & \ 4 \epsilon \left( \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right) - V_{LJ}(r_c) - (r - r_c) \frac{\partial}{\partial r} (V_{LJ}(r)) \Big\vert_{r_c} \ \ \ \  if \  r  \leqslant r_c \\
	V_{LJ}^{ts}(r) \ & \ 0      \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \;\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; if \  r  > r_c \\
\end{array}
\right.
\]
$\\ \\$
It is the same potential that Lennard-Jonnes, but we consider that the particle which are at a distance superior at $r_c$ don't influence the force. This reduce a lot the computer calcul time especially for lots of particles.

\section{Algorithms}

Now it needs one algorithm to calculate position, velocity and acceleration of each particle at each step of time (as we consider a disceet system), so we will used here two different algorithms which is Velocity Verlet algorithm and Beeman algorithm. Don't forget that acceleration is deduce by the force, and the force is calculated by the potential, and the potential is given by the position (and the potential kind considered).

\subsection{Velocity Verlet}

1. Calculate:
\begin{equation}
	 \vec{v}(t + \frac{1}{2}\Delta t) = \vec{v}(t) + \frac{1}{2}\vec{a}(t)\Delta t
\end{equation}

2. Calculate:
\begin{equation}
 	 \vec{x}(t + \Delta t) = \vec{x}(t) + \vec{v} (t + \frac{1}{2} \Delta t) \Delta t
\end{equation}
 \\

3. Derive $\vec{a}(t + \Delta t)$ from the interaction potential using $\vec{x}(t + \Delta t)$
\\ 

4. Calculate:
\begin{equation}
	 \vec{v}(t + \Delta t) = \vec{v}(t + \frac{1}{2}\Delta t) + \frac{1}{2}\vec{a}(t + \Delta t)\Delta t
\end{equation}
\\
	
The velocity at the following step of time is now calculated.

\subsection{Beeman}
 
\begin{equation}
    \vec{x}(t+\Delta t) = \vec{x}(t) + \vec{v}(t) \Delta t + \left(\frac{2}{3}\vec{a}(t) - \frac{1}{6} \vec{a}(t - \Delta t) \right)\Delta t^2 
\end{equation}

\begin{equation}
    \vec{v}(t + \Delta t) = \vec{v}(t) + \left(\frac{1}{3}\vec{a}(t + \Delta t) + \frac{5}{6}\vec{a}(t) - \frac{1}{6}\vec{a}(t - \Delta t) \right) \Delta t
\end{equation}
\\

Those equation depends on the acceleration at the precedent step of time.
This algorithm is at the third order, so with more precision than Velocity Verlet (second order). 



\section{Temperature and Pressure}

The temperature is called kinetic temperature because implied by particles movments, so we can deduce kinetic temperature of the kinetic energy(where N is atom number, and $k_b$ the Boltzmann constant)  %thermalisation
\begin{equation}
	\boxed{T = \frac{2}{3 N k_b}  \sum ^N_{i=1} m_i v_i^2}
\end{equation}

The pressure will calculate in the part of the potential because has two composants, one kinetic and the second called pressure of excess (virial equation), this last is emplied by the force and the position between particles.
\begin{equation}
	\boxed{P = \frac{N}{V}  k_b  T + \frac{1}{3V}  \sum ^N_{i=1} \sum ^N_{j>i} f(r_{ij}) . r_{ij}}
\end{equation}

\section{Simulation Units}

In general when make simulation, we search to express each variable as a variable without dimension. So we want write each variable as a function of $\sigma$ (distance at which the inter-particle potential is zero), $\epsilon$ (depth of the potential well) and $m$ (mass of one atom) which will rescpectively in our program the units of lenght, energy and mass. 
We can deduce without dimension variable:


-Energy: \begin{equation}
		    \tilde{E} = \frac{E}{\epsilon} 
		 \end{equation}

-Volume: \begin{equation}
		 	\tilde{V} = \frac{V}{\sigma^3} 
		 \end{equation}
		 
-Time:   \begin{equation}
			\tilde{t} = \frac{t}{\sigma} \sqrt{\frac{\epsilon}{m}}
		 \end{equation}

-Temperature: \begin{equation}
		 	  	\tilde{T} = \frac{k_b T}{\epsilon} 
		      \end{equation}

-Pressure: \begin{equation}
		   		\tilde{P} = \frac{\sigma^3}{\epsilon} P
		   \end{equation}

-Density: \begin{equation}
		 	\tilde{\rho} = \sigma^3 \rho 
		  \end{equation}

\section{Phase diagram}
We show a qualitative phase diagram:

%\begin{figure}[htbp] %  figure placement: here, top, bottom, or page
\begin{center}
	   \includegraphics[width=0.5 \columnwidth]{phaselj.png} 
	   \label{fig:figure5}
\end{center}
%\end{figure}

To Lennard-Jones potential, by Molecular Dynamics and Monte-Carlo simulations we find $\tilde{T}_c = 1.313$, $\tilde{\rho_c} = 0.310$ and $\tilde{T_t} = 0.69$ (where $\tilde{T}_c$, $\tilde{\rho_c}$ and $\tilde{T_t}$ are in reduced units).

To truncated and shifted Lennard-Jonnes potential, we find $\tilde{T}_c = 1.085$ and $\tilde{\rho_c} = 0.317$ (for $\sigma_c = 2.5$, where $\sigma_c$ is the beam of truncation).

\section{CPU time and optimization methods}
A point which is very important in molecular dynamic simulation is the "CPU time", because more atoms there are, the more computer calcul time is long. In fact that take the bigest time at the computer into the program is the double loop over atoms, so  during the force and potential calculation, we can express the calculation time simply by a qualitative method by couting loop number as function of atom number.

\subsubsection{General method}
In general, as we discuss previously, we have just one box (simulation box which have a lenght L) and we must take each atoms and calculate there interaction potential with all the others, as it's shown on the following figure.

$\\ \\$

%\begin{figure}[htbp] %  figure placement: here, top, bottom, or page
\begin{center}
	   \includegraphics[width=0.5 \columnwidth]{resultat3.pdf}
	   \label{fig:figure8}
\end{center}
%\end{figure}

$\\ \\ \\$

This method during: 

\begin{equation}
	\frac{N(N-1)}{2}
	\label{eq1}
\end{equation}

Thus the time is $N^2$, and if atom number is very important like one million or more (it's frequent in molecular dynamics) the time is very very long, because depends on atoms number square.

\subsubsection{Verlet cells}
We must use truncated potential to this method. This method allows of not test all the others particles but only those are at a distance approximately equal at $r_c$. We take generaly $\sigma_c=2.7$, where $\sigma_c$ is the lenght of one cell. We must cut simulation box at a distance approximately equal at 2.7 and test only neighbour cells, that reduce a lot calculation time. Don't forget the PBC. We show Verlet cells on the following figure.

$\\ \\ \\$

%\begin{figure}[htbp] %  figure placement: here, top, bottom, or page
\begin{center}
	   \includegraphics[width=0.5 \columnwidth]{resultat4.pdf}
	   \label{fig:figure8}
\end{center}
%\end{figure}

$\\ \\ \\$

We consider in general M cells, that imply $\frac{N}{M}$ atoms by cells on average.
This method during:

\begin{equation}
	 M \left(\frac{\frac{N}{M}\left(\frac{N}{M}-1\right)}{2} + 26 \left(\frac{N}{M}\right)^2\right)
\end{equation}

The first therm correspond at the interaction between particles which is in the same cell and the second therm correspond at the interaction between the particles of one cell with the particle of another, multiplied by M because calculate for every cells. We find after development:

\begin{equation}
	 \frac{1}{2}\left(\frac{53 N^2}{M} - N\right)
	 \label{eq2}
\end{equation}

And we know the following relation:  $N_c = \frac{L}{\sigma_c}$, $M = N_c^3=\left(\frac{L}{\sigma_c}\right)^3$, $\rho = \frac{N}{L^3}$ (where $N_c$ is cells number by dimension).

We can deduce another expression for (\ref{eq2}) another expression as function of density and atoms number:

\begin{equation}
	 \frac{N}{2} \left(53 \rho \sigma_c^3 - 1\right)
	 \label{eq3}
\end{equation}

We shall set afterward to illustrate that curve a value for $\sigma_c$ which is the same that given before. So: $53 \sigma_c^3 \simeq 830 $.

$\\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\$

So we compare the two methods in ploting on graphics:

%\begin{figure}[htbp] %  figure placement: here, top, bottom, or page
\begin{center}
	   \includegraphics[width=0.5 \columnwidth]{resultat.png}
	   \label{fig:figure6}
\end{center}
%\end{figure}

%\begin{figure}[htbp] %  figure placement: here, top, bottom, or page
\begin{center}
	   \includegraphics[width=0.5 \columnwidth]{resultat2.png} 
	   \label{fig:figure7}
\end{center}
%\end{figure}

We see to a same atoms number $10^6$, for a solid (which have the bigest density $\simeq$ 30) too, Verlet cells methods make one thousand time less loops than general method.


\end{document}  